home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / derfc.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  12.7 KB  |  212 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((nterf 0)
  12.       (nterfc 0)
  13.       (nterc2 0)
  14.       (xsml 0.0)
  15.       (xmax 0.0)
  16.       (sqeps 0.0)
  17.       (erfcs (make-array 21 :element-type 'double-float))
  18.       (erc2cs (make-array 49 :element-type 'double-float))
  19.       (erfccs (make-array 59 :element-type 'double-float))
  20.       (sqrtpi 1.772453850905516)
  21.       (first nil))
  22.   (declare (type f2cl-lib:logical first)
  23.            (type (simple-array double-float (59)) erfccs)
  24.            (type (simple-array double-float (49)) erc2cs)
  25.            (type (simple-array double-float (21)) erfcs)
  26.            (type double-float sqrtpi sqeps xmax xsml)
  27.            (type f2cl-lib:integer4 nterc2 nterfc nterf))
  28.   (f2cl-lib:fset (f2cl-lib:fref erfcs (1) ((1 21))) -0.04904612123469181)
  29.   (f2cl-lib:fset (f2cl-lib:fref erfcs (2) ((1 21))) -0.14226120510371365)
  30.   (f2cl-lib:fset (f2cl-lib:fref erfcs (3) ((1 21))) 0.010035582187599796)
  31.   (f2cl-lib:fset (f2cl-lib:fref erfcs (4) ((1 21))) -5.768764699767486e-4)
  32.   (f2cl-lib:fset (f2cl-lib:fref erfcs (5) ((1 21))) 2.7419931252196067e-5)
  33.   (f2cl-lib:fset (f2cl-lib:fref erfcs (6) ((1 21))) -1.1043175507344509e-6)
  34.   (f2cl-lib:fset (f2cl-lib:fref erfcs (7) ((1 21))) 3.8488755420345033e-8)
  35.   (f2cl-lib:fset (f2cl-lib:fref erfcs (8) ((1 21))) -1.1808582533875464e-9)
  36.   (f2cl-lib:fset (f2cl-lib:fref erfcs (9) ((1 21))) 3.233421582605091e-11)
  37.   (f2cl-lib:fset (f2cl-lib:fref erfcs (10) ((1 21))) -7.991015947004547e-13)
  38.   (f2cl-lib:fset (f2cl-lib:fref erfcs (11) ((1 21))) 1.7990725113961456e-14)
  39.   (f2cl-lib:fset (f2cl-lib:fref erfcs (12) ((1 21))) -3.7186354878186934e-16)
  40.   (f2cl-lib:fset (f2cl-lib:fref erfcs (13) ((1 21))) 7.103599003714253e-18)
  41.   (f2cl-lib:fset (f2cl-lib:fref erfcs (14) ((1 21))) -1.2612455119155225e-19)
  42.   (f2cl-lib:fset (f2cl-lib:fref erfcs (15) ((1 21))) 2.0916406941769294e-21)
  43.   (f2cl-lib:fset (f2cl-lib:fref erfcs (16) ((1 21))) -3.2539731029314073e-23)
  44.   (f2cl-lib:fset (f2cl-lib:fref erfcs (17) ((1 21))) 4.7668672097976744e-25)
  45.   (f2cl-lib:fset (f2cl-lib:fref erfcs (18) ((1 21))) -6.598012078285136e-27)
  46.   (f2cl-lib:fset (f2cl-lib:fref erfcs (19) ((1 21))) 8.655011469963763e-29)
  47.   (f2cl-lib:fset (f2cl-lib:fref erfcs (20) ((1 21))) -1.0788925177498063e-30)
  48.   (f2cl-lib:fset (f2cl-lib:fref erfcs (21) ((1 21))) 1.2811883993017004e-32)
  49.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (1) ((1 49))) -0.0696013466023095)
  50.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (2) ((1 49))) -0.04110133936262089)
  51.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (3) ((1 49))) 0.003914495866689627)
  52.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (4) ((1 49))) -4.906395650548979e-4)
  53.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (5) ((1 49))) 7.157479001377036e-5)
  54.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (6) ((1 49))) -1.153071634131233e-5)
  55.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (7) ((1 49))) 1.994670590201998e-6)
  56.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (8) ((1 49))) -3.642666471599223e-7)
  57.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (9) ((1 49))) 6.944372610005013e-8)
  58.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (10) ((1 49))) -1.3712209021043661e-8)
  59.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (11) ((1 49))) 2.7883896610071374e-9)
  60.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (12) ((1 49))) -5.814164724331161e-10)
  61.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (13) ((1 49))) 1.2389204917527533e-10)
  62.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (14) ((1 49))) -2.690639145306744e-11)
  63.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (15) ((1 49))) 5.94261435084791e-12)
  64.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (16) ((1 49))) -1.3323867357581193e-12)
  65.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (17) ((1 49))) 3.0280468061771326e-13)
  66.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (18) ((1 49))) -6.966648814941033e-14)
  67.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (19) ((1 49))) 1.6208545410539232e-14)
  68.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (20) ((1 49))) -3.809934465250492e-15)
  69.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (21) ((1 49))) 9.040487815978833e-16)
  70.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (22) ((1 49))) -2.1640061950896078e-16)
  71.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (23) ((1 49))) 5.222102233995854e-17)
  72.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (24) ((1 49))) -1.2697296023645555e-17)
  73.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (25) ((1 49))) 3.109145504276198e-18)
  74.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (26) ((1 49))) -7.663762920320385e-19)
  75.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (27) ((1 49))) 1.9008192513627456e-19)
  76.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (28) ((1 49))) -4.7422072790690395e-20)
  77.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (29) ((1 49))) 1.1896492000765284e-20)
  78.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (30) ((1 49))) -3.00003559032578e-21)
  79.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (31) ((1 49))) 7.602993453043245e-22)
  80.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (32) ((1 49))) -1.9359094476068728e-22)
  81.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (33) ((1 49))) 4.9513991247733385e-23)
  82.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (34) ((1 49))) -1.271807481336372e-23)
  83.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (35) ((1 49))) 3.280049600469513e-24)
  84.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (36) ((1 49))) -8.492320176822895e-25)
  85.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (37) ((1 49))) 2.20691789280756e-25)
  86.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (38) ((1 49))) -5.755617245696529e-26)
  87.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (39) ((1 49))) 1.5061915336392342e-26)
  88.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (40) ((1 49))) -3.954502959018798e-27)
  89.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (41) ((1 49))) 1.0415297041515008e-27)
  90.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (42) ((1 49))) -2.751487795278766e-28)
  91.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (43) ((1 49))) 7.290058205497557e-29)
  92.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (44) ((1 49))) -1.9369396459159477e-29)
  93.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (45) ((1 49))) 5.160357112051487e-30)
  94.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (46) ((1 49))) -1.3784193221930938e-30)
  95.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (47) ((1 49))) 3.691326793107069e-31)
  96.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (48) ((1 49))) -9.909389590624368e-32)
  97.   (f2cl-lib:fset (f2cl-lib:fref erc2cs (49) ((1 49))) 2.666491705195389e-32)
  98.   (f2cl-lib:fset (f2cl-lib:fref erfccs (1) ((1 59))) 0.07151793102029248)
  99.   (f2cl-lib:fset (f2cl-lib:fref erfccs (2) ((1 59))) -0.026532434337606714)
  100.   (f2cl-lib:fset (f2cl-lib:fref erfccs (3) ((1 59))) 0.001711153977920856)
  101.   (f2cl-lib:fset (f2cl-lib:fref erfccs (4) ((1 59))) -1.6375166345851788e-4)
  102.   (f2cl-lib:fset (f2cl-lib:fref erfccs (5) ((1 59))) 1.987129350055204e-5)
  103.   (f2cl-lib:fset (f2cl-lib:fref erfccs (6) ((1 59))) -2.843712412766556e-6)
  104.   (f2cl-lib:fset (f2cl-lib:fref erfccs (7) ((1 59))) 4.60616130896313e-7)
  105.   (f2cl-lib:fset (f2cl-lib:fref erfccs (8) ((1 59))) -8.227753025879208e-8)
  106.   (f2cl-lib:fset (f2cl-lib:fref erfccs (9) ((1 59))) 1.5921418727709008e-8)
  107.   (f2cl-lib:fset (f2cl-lib:fref erfccs (10) ((1 59))) -3.295071362252843e-9)
  108.   (f2cl-lib:fset (f2cl-lib:fref erfccs (11) ((1 59))) 7.223439760400555e-10)
  109.   (f2cl-lib:fset (f2cl-lib:fref erfccs (12) ((1 59))) -1.6648558133987296e-10)
  110.   (f2cl-lib:fset (f2cl-lib:fref erfccs (13) ((1 59))) 4.010392588237665e-11)
  111.   (f2cl-lib:fset (f2cl-lib:fref erfccs (14) ((1 59))) -1.0048162144257311e-11)
  112.   (f2cl-lib:fset (f2cl-lib:fref erfccs (15) ((1 59))) 2.6082759133003336e-12)
  113.   (f2cl-lib:fset (f2cl-lib:fref erfccs (16) ((1 59))) -6.991110560404025e-13)
  114.   (f2cl-lib:fset (f2cl-lib:fref erfccs (17) ((1 59))) 1.929492333261707e-13)
  115.   (f2cl-lib:fset (f2cl-lib:fref erfccs (18) ((1 59))) -5.470131188754331e-14)
  116.   (f2cl-lib:fset (f2cl-lib:fref erfccs (19) ((1 59))) 1.5896633097626972e-14)
  117.   (f2cl-lib:fset (f2cl-lib:fref erfccs (20) ((1 59))) -4.726893980197554e-15)
  118.   (f2cl-lib:fset (f2cl-lib:fref erfccs (21) ((1 59))) 1.435873376784985e-15)
  119.   (f2cl-lib:fset (f2cl-lib:fref erfccs (22) ((1 59))) -4.4495105618173586e-16)
  120.   (f2cl-lib:fset (f2cl-lib:fref erfccs (23) ((1 59))) 1.4048108847682333e-16)
  121.   (f2cl-lib:fset (f2cl-lib:fref erfccs (24) ((1 59))) -4.51381838776421e-17)
  122.   (f2cl-lib:fset (f2cl-lib:fref erfccs (25) ((1 59))) 1.4745215410451332e-17)
  123.   (f2cl-lib:fset (f2cl-lib:fref erfccs (26) ((1 59))) -4.8926214069457763e-18)
  124.   (f2cl-lib:fset (f2cl-lib:fref erfccs (27) ((1 59))) 1.647612141410647e-18)
  125.   (f2cl-lib:fset (f2cl-lib:fref erfccs (28) ((1 59))) -5.62681717632941e-19)
  126.   (f2cl-lib:fset (f2cl-lib:fref erfccs (29) ((1 59))) 1.9474433822320786e-19)
  127.   (f2cl-lib:fset (f2cl-lib:fref erfccs (30) ((1 59))) -6.82630564294842e-20)
  128.   (f2cl-lib:fset (f2cl-lib:fref erfccs (31) ((1 59))) 2.4219888872986495e-20)
  129.   (f2cl-lib:fset (f2cl-lib:fref erfccs (32) ((1 59))) -8.693414133503069e-21)
  130.   (f2cl-lib:fset (f2cl-lib:fref erfccs (33) ((1 59))) 3.1551803462280853e-21)
  131.   (f2cl-lib:fset (f2cl-lib:fref erfccs (34) ((1 59))) -1.1573723240496088e-21)
  132.   (f2cl-lib:fset (f2cl-lib:fref erfccs (35) ((1 59))) 4.288947161605653e-22)
  133.   (f2cl-lib:fset (f2cl-lib:fref erfccs (36) ((1 59))) -1.6050307420576168e-22)
  134.   (f2cl-lib:fset (f2cl-lib:fref erfccs (37) ((1 59))) 6.063298757453802e-23)
  135.   (f2cl-lib:fset (f2cl-lib:fref erfccs (38) ((1 59))) -2.3114042516979588e-23)
  136.   (f2cl-lib:fset (f2cl-lib:fref erfccs (39) ((1 59))) 8.888778540661886e-24)
  137.   (f2cl-lib:fset (f2cl-lib:fref erfccs (40) ((1 59))) -3.4472605766513764e-24)
  138.   (f2cl-lib:fset (f2cl-lib:fref erfccs (41) ((1 59))) 1.347865460206965e-24)
  139.   (f2cl-lib:fset (f2cl-lib:fref erfccs (42) ((1 59))) -5.311794071125021e-25)
  140.   (f2cl-lib:fset (f2cl-lib:fref erfccs (43) ((1 59))) 2.1093410586197825e-25)
  141.   (f2cl-lib:fset (f2cl-lib:fref erfccs (44) ((1 59))) -8.43836558792379e-26)
  142.   (f2cl-lib:fset (f2cl-lib:fref erfccs (45) ((1 59))) 3.3999825249452087e-26)
  143.   (f2cl-lib:fset (f2cl-lib:fref erfccs (46) ((1 59))) -1.3794523880732423e-26)
  144.   (f2cl-lib:fset (f2cl-lib:fref erfccs (47) ((1 59))) 5.634490311833252e-27)
  145.   (f2cl-lib:fset (f2cl-lib:fref erfccs (48) ((1 59))) -2.3164904344770657e-27)
  146.   (f2cl-lib:fset (f2cl-lib:fref erfccs (49) ((1 59))) 9.58446284460181e-28)
  147.   (f2cl-lib:fset (f2cl-lib:fref erfccs (50) ((1 59))) -3.9907228803301104e-28)
  148.   (f2cl-lib:fset (f2cl-lib:fref erfccs (51) ((1 59))) 1.6721292259444773e-28)
  149.   (f2cl-lib:fset (f2cl-lib:fref erfccs (52) ((1 59))) -7.045991522766013e-29)
  150.   (f2cl-lib:fset (f2cl-lib:fref erfccs (53) ((1 59))) 2.9797684028642063e-29)
  151.   (f2cl-lib:fset (f2cl-lib:fref erfccs (54) ((1 59))) -1.262522466460619e-29)
  152.   (f2cl-lib:fset (f2cl-lib:fref erfccs (55) ((1 59))) 5.395438704542488e-30)
  153.   (f2cl-lib:fset (f2cl-lib:fref erfccs (56) ((1 59))) -2.380992882531459e-30)
  154.   (f2cl-lib:fset (f2cl-lib:fref erfccs (57) ((1 59))) 1.0990528301027616e-30)
  155.   (f2cl-lib:fset (f2cl-lib:fref erfccs (58) ((1 59))) -4.867713741644965e-31)
  156.   (f2cl-lib:fset (f2cl-lib:fref erfccs (59) ((1 59))) 1.5258772641103577e-31)
  157.   (setq first f2cl-lib:%true%)
  158.   (defun derfc (x)
  159.     (declare (type double-float x))
  160.     (prog ((txmax 0.0) (y 0.0) (derfc 0.0) (eta 0.0f0))
  161.       (declare (type single-float eta) (type double-float derfc y txmax))
  162.       (cond
  163.        (first (setf eta (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3))))
  164.               (setf nterf (initds erfcs 21 eta))
  165.               (setf nterfc (initds erfccs 59 eta))
  166.               (setf nterc2 (initds erc2cs 49 eta))
  167.               (setf xsml
  168.                       (-
  169.                        (f2cl-lib:fsqrt
  170.                         (- (f2cl-lib:flog (* sqrtpi (f2cl-lib:d1mach 3)))))))
  171.               (setf txmax
  172.                       (f2cl-lib:fsqrt
  173.                        (- (f2cl-lib:flog (* sqrtpi (f2cl-lib:d1mach 1))))))
  174.               (setf xmax
  175.                       (- (+ txmax (/ (* -0.5 (f2cl-lib:flog txmax)) txmax))
  176.                          0.01))
  177.               (setf sqeps (f2cl-lib:fsqrt (* 2.0 (f2cl-lib:d1mach 3))))))
  178.       (setf first f2cl-lib:%false%)
  179.       (if (> x xsml) (go label20))
  180.       (setf derfc 2.0)
  181.       (go end_label)
  182.      label20
  183.       (if (> x xmax) (go label40))
  184.       (setf y (coerce (abs x) 'double-float))
  185.       (if (> y 1.0) (go label30))
  186.       (if (< y sqeps) (setf derfc (+ 1.0 (/ (* -2.0 x) sqrtpi))))
  187.       (if (>= y sqeps)
  188.           (setf derfc
  189.                   (- 1.0
  190.                      (* x (+ 1.0 (dcsevl (- (* 2.0 x x) 1.0) erfcs nterf))))))
  191.       (go end_label)
  192.      label30
  193.       (setf y (* y y))
  194.       (if (<= y 4.0)
  195.           (setf derfc
  196.                   (* (/ (exp (- y)) (abs x))
  197.                      (+ 0.5
  198.                         (dcsevl (/ (- (/ 8.0 y) 5.0) 3.0) erc2cs nterc2)))))
  199.       (if (> y 4.0)
  200.           (setf derfc
  201.                   (* (/ (exp (- y)) (abs x))
  202.                      (+ 0.5 (dcsevl (- (/ 8.0 y) 1.0) erfccs nterfc)))))
  203.       (if (< x 0.0) (setf derfc (- 2.0 derfc)))
  204.       (go end_label)
  205.      label40
  206.       (xermsg "SLATEC" "DERFC" "X SO BIG ERFC UNDERFLOWS" 1 1)
  207.       (setf derfc 0.0)
  208.       (go end_label)
  209.      end_label
  210.       (return (values derfc nil)))))
  211.  
  212.